This is what ought to be the reproduction of the analyses presented in the article Inferring biotic interactions from proxies by Ignacio Morales-Castilla, Miguel G. Matias, Dominique Gravel, and Miguel B. Ara??jo (2015, TREE). The article can be found here. Unfortuntely, reproducing thier analyses requires to go deep into two other studies, one about the Serengeti food web (Baskerville et al. (2011, PLoS Computational Biology)), and one about three Caribean coral reefs (Roopnarine & Hertog (2013, Dataset Papers in Ecology)).
Clear R environment and load libraries.
# Clear workspace.
rm(list=ls())
# Load libraries.
library(igraph)
##
## Attaching package: 'igraph'
##
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
##
## The following object is masked from 'package:base':
##
## union
library(RCurl)
## Loading required package: bitops
library(readxl)
Load the data from GitHub straight into R. This is more convenient for everyone because one does not have to always change the local file path…
# Species in the Serengeti food web.
dd <- getURL("https://raw.githubusercontent.com/opetchey/RREEBES/master/MORALES-CASTILLA_etal_2015_TREE/Data/Serengeti%20data/journal.pcbi.1002321.s004.CSV",
cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
dd.01 <- read.csv(text = dd)
# Feeding links in the Serengeti food web.
dd <- getURL("https://raw.githubusercontent.com/opetchey/RREEBES/master/MORALES-CASTILLA_etal_2015_TREE/Data/Serengeti%20data/journal.pcbi.1002321.s005.CSV",
cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
dd.02 <- read.csv(text = dd)
# Consensus partition.
dd <- getURL("https://raw.githubusercontent.com/opetchey/RREEBES/master/MORALES-CASTILLA_etal_2015_TREE/Data/Serengeti%20data/journal.pcbi.1002321.s006.CSV",
cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
dd.03 <- read.csv(text = dd)
# Link densities between groups in the consensus partition.
dd <- getURL("https://raw.githubusercontent.com/opetchey/RREEBES/master/MORALES-CASTILLA_etal_2015_TREE/Data/Serengeti%20data/journal.pcbi.1002321.s007.CSV",
cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
dd.04 <- read.csv(text = dd)
Calculate the number of species and the number of potential links one could expect from this.
# Species only in Cuban food web.
n.species.Serengeti <- length(dd.01[,1])
n.links.species.Serengeti <- n.species.Serengeti^2-n.species.Serengeti
Note: There are 175 species in the Serengeti food web. This results in 25760 potential links (\(n^2-n\)).
dd <- getURL("https://raw.githubusercontent.com/opetchey/RREEBES/Morales-Castilla_etal_2015_TREE/MORALES-CASTILLA_etal_2015_TREE/Data/Coral%20reef%20data/857470.item.4%20-%20Guild%20key%2C%20all%20data%20sets.csv",
cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
dd.05 <- read.csv(text = dd)
str(dd.05)
## 'data.frame': 755 obs. of 6 variables:
## $ Guild.Number : int 1 2 2 2 2 2 2 2 2 2 ...
## $ Taxa : Factor w/ 753 levels ""," ",..: 1 661 120 126 127 132 282 396 407 685 ...
## $ Fish.Body.Length: num NA NA NA NA NA NA NA NA NA NA ...
## $ Cayman.Islands : Factor w/ 3 levels ""," ","x ": 1 1 1 1 1 1 1 1 1 1 ...
## $ Cuba : Factor w/ 3 levels ""," ","x": 1 1 1 1 1 1 1 1 1 1 ...
## $ Jamaica : Factor w/ 3 levels ""," ","x": 1 1 1 1 1 1 1 1 1 1 ...
# Only data with species from the Cuban food web.
# Guilds only in Cuban food web.
dd.05.a <- subset(dd.05, Cuba=="x")
# Only the fish pecies in Cuban food web.
n.species.fish.Cuba <- length(subset(dd.05, Cuba=="x")$Cuba)
# Trophic data for the Cuban food web
dd <- getURL("https://raw.githubusercontent.com/opetchey/RREEBES/Morales-Castilla_etal_2015_TREE/MORALES-CASTILLA_etal_2015_TREE/Data/Coral%20reef%20data/857470.item.2%20-%20Trophic%20data%20for%20Cuba_sorted.csv",
cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
dd.06 <- read.csv(text = dd)
# Guilds in Cuban food web.
n.guilds.Cuba <- length(dd.06[,1])
n.links.guilds.Cuba <- n.guilds.Cuba^2-n.guilds.Cuba
“The procedure for building the backbone of an interaction network starts with the identification of species more likely to share similar interactions. The concept is similar to that of modules [28], but we avoid this terminology because modules are usually determined a posteriori and can also refer to simple assemblages of species such as linear food chains or apparent competition [29]. Instead, we define interacting groups based on a priori expectations of interactions. The concept is also analogous to that of guilds [30]. Guilds, however, are restricted to species shar- ing similar resources, and thus do not encompass non- consumptive interactions such as competition or niche construction [27]. A flexible definition of interacting groups based on traits, phylogenies, and geographical distribu- tions would enable combination of heterogeneous information.” (page 4, bottom-right)
“Interacting groups of species are defined a priori to simplify the removal of forbidden links. The groups were defined based on the trophic hierarchy of the different species within each eco- system (e.g., primary producers, grazers, small and large carnivores). This process of trophic classification of species led to identification of forbidden links and removal of 30% of all potential direct links in the coral reef, and 22% in the Serengeti (e.g., herbivores eating predators; Figure 1).” (page 5, bottom-left)
“Refinement of the species groupings was achieved by con- sidering the characteristics of the consumer species (e.g., distinguishing small versus large carnivores in the Seren- geti example, or separating invertebrate feeders, omnivo- rous, and carnivorous fish in the coral reef example; Figure 1).” (page 5, bottom-right)
Taking all potential links into account the food webs would look like this:
dd.01.a <- expand.grid(dd.01$code, dd.01$code)
# Build the graph.
f <- graph.data.frame(dd.01.a, directed=F, vertices=NULL)
# Now, remove inraspecific trophic interactions - there is no cannibalism.
dd.01.b <- dd.01.a[which(dd.01.a[,1]!=dd.01.a[,2]),]
# Now, plot the graph.
plot.igraph(f,
layout=layout.circle,
main="Serengeti food web (no information)",
vertex.size=8,
vertex.color=adjustcolor("#A6B07E", alpha.f=0.6),
vertex.label.cex=0.5,
vertex.label.family = "sans",
edge.width=0.2,
edge.color=adjustcolor("#999999", alpha.f=0.3),
vertex.label.color="#7D1D3F",
vertex.frame.color="white")
This is what the real food web looks like. So far, there is no ordering according to any traits (e.g., nutrition, size…)
# Build the graph.
g <- graph.data.frame(dd.02, directed=F, vertices=NULL)
# Display the food web graph.
# tkplot(g)
# plot.igraph(g)
plot.igraph(g,
# layout=layout.fruchterman.reingold,
layout=layout.circle,
main="Serengeti food web (real food web)",
vertex.size=8,
vertex.color=adjustcolor("#A6B07E", alpha.f=0.6),
vertex.label.cex=0.5,
vertex.label.family = "sans",
edge.width=0.2,
edge.color=adjustcolor("#333333", alpha.f=0.6),
vertex.label.color="#7D1D3F",
vertex.frame.color="white")
Once we know about the trophic levels the species/guilds belong to we can use layout.reingold.tilford(g, flip.y=FALSE) to introduce some hierarchy. For example, if we assign three trophic levels just randomly then the food web can be displayed like this:
# plot.igraph(g,
# # layout=layout.fruchterman.reingold,
# layout=layout.reingold.tilford(g, flip.y=FALSE),
# main="Serengeti food web (no trait information)",
# vertex.size=8,
# vertex.color=adjustcolor("#A6B07E", alpha.f=0.6),
# vertex.label.cex=0.5,
# vertex.label.family = "sans",
# edge.width=2,
# vertex.label.color="#7D1D3F",
# vertex.frame.color="white")
dd.06.a <- expand.grid(dd.06$Guild.Description, dd.06$Guild.Description)
# Build the graph.
f <- graph.data.frame(dd.06.a, directed=F, vertices=NULL)
# Now, remove inraspecific trophic interactions - there is no cannibalism.
dd.06.b <- dd.06.a[which(dd.06.a[,1]!=dd.06.a[,2]),]
# Should equal the number of potential links, 53.130.
length(dd.06.b[,1])
## [1] 53128
# Now, plot the graph.
plot.igraph(f,
layout=layout.circle,
main="Cuban coral reef food web (no information)",
vertex.size=8,
vertex.color=adjustcolor("#A6B07E", alpha.f=0.6),
vertex.label.cex=0.5,
vertex.label.family = "sans",
edge.width=0.2,
edge.color=adjustcolor("#999999", alpha.f=0.3),
vertex.label.color="#7D1D3F",
vertex.frame.color="white")
# If number of prey equals 0 than the species is a primary producer.
Comments, errors & modifications
Article text
Data sets